Computational complexity of quantum optimal control landscapes 
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We study the Hamiltonian-independent contribution to the complexity of quantum optimal con- 
trol problems. The optimization of controls that steer quantum systems to desired objectives can 
itself be considered a classical dynamical system that executes an analog computation. The system- 
independent component of the equations of motion of this dynamical system can be integrated 
analytically for various classes of discrete quantum control problems. For the maximization of ob- 
servable expectation values from an initial pure state and the maximization of the fidelity of quantum 
gates, the time complexity of the corresponding computation belongs to the class continuous log 
(CLOG), the lowest analog complexity class, equivalent to the discrete complexity class NC. The 
simple scaling of the Hamiltonian-independent contribution to these problems with quantum system 
dimension indicates that with appropriately designed search algorithms, quantum optimal control 
can be rendered efficient even for large systems. 
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I. INTRODUCTION 

In recent years the methodology of optimal control 
theory has been applied to achieve various objectives in 
quantum systems. The two most common objectives are 
the maximization of the expectation value of an observ- 
able, and the maximization of the fidelity of a unitary 
transformation of the wavefunction. Whereas the maxi- 
mization of observable expectation values has met with 
widespread success in experimental and computational 
incarnations [l[ , the achievement of high fidelity unitary 
transformations has proven more challenging [2|], espe- 
cially for large systems. However, it is not clear whether 
these apparent distinctions are specific to the algorithms 
employed, or whether they are indicative of universal 
and inherent features of the quantum optimization prob- 
lems. 

Given the computational expense of propagating the 
Schrodinger equation, the scaling with system size of 
the expense of executing a search over the space of con- 
trol fields - referred to as the problem's computational 
complexity - plays an essential role in determining the 
feasibility and efficiency of control optimizations in large 
quantum systems. The complexity of control optimiza- 
tion for quantum gates is of particular interest, as it 
corresponds to the classical complexity associated with 
the physical implementation of quantum logical oper- 
ations. Recent studies have addressed the computa- 
tional complexity (and related, the critical topology) 
of several associated problems in quantum information 
science, including quantum state and process identifi- 
cation In contrast to these tasks, the variational 
problem of quantum control optimization occurs on an 
infinite-dimensional parameter space. However, for dis- 
crete quantum systems, both the complexity and critical 
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topology of quantum control problems can be divided 
into finite- and infinite-dimensional contributions, sim- 
plifying their study. 

Consider first the critical topology of quantum control 
variational problems. Denoting the cost functional by J 
and the control field by e(t), we have according to the 
chain rule 

SJ dJ SU(T) 

~ dU(T) ' Se(t) ' 

where U(T) is the dynamical propagator resulting from 
application of the control e(t) over the time interval 
< t < T. As such, the critical points of quantum con- 
trol problems fall into two classes. The first type of min- 
imizer corresponds to those control Hamiltonians that 
are critical points of the control objective functional, 
but are not critical points of the map between control 
fields and associated dynamical propagators (i.e., points 
at which d ^'^ = 0, while the Frechet derivative map- 
ping from the control variation Se(t) to SU(T) at t = T 
is surjective). These critical points are called normal 
extremal controls. The second type of minimizer corre- 
sponds to critical points of the latter map (i.e., points 
at which the mapping from Se(t) to SU(T) is not locally 
surjective). These latter points are referred to as ab- 
normal extremal controls. For several low-dimensional 
quantum control problems, it has been shown that ab- 
normal extremals are particularly scarce 0, H, @| , sug- 
gesting that the normal extremal controls dominate the 
critical point topology of these problems. 

As such, attention has focused on characterizing the 
system-independent (or kinematic) contribution to the 
critical topology of quantum optimal control problems 
@j B ■ For a given quantum system, several different 
classes of quantum control problems may be envisioned, 
each of which has its own associated kinematic critical 
topology. Besides the distinction between quantum gate 
and observable control, observable control itself may be 
subdivided into several qualitatively distinct optimiza- 
tion problems, depending on the initial state (density 
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matrix) of the system and the type of observable to be 
maximized. Mathematically, the critical topologies of 
observable control problems can be classified on the ba- 
sis of the eigenvalue spectra of the initial density matrix 
po and the observable operator <d @. To date, however, 
the system-independent contribution to the computa- 
tional complexity of these problems has not been ad- 
dressed. Here, we show that for several classes of quan- 
tum control problems, this complexity can, remarkably, 
be analytically determined. 



II. CONTROL OPTIMIZATION AS A 
DYNAMICAL SYSTEM 

A universal quantum optimal control cost functional 
can be written as: 



J = $(J7(T),T) 
ReTr ' 



{^ + \H{eit))U(t))Bu s \,U 



\f^\e{t)\ 2 dt (1) 
Jo s{t) 



where B(t) is a Lagrange multiplier constraining the 
quantum system dynamics to obey the Schrodinger 
equation, e(t) denotes the time-dependent control field, 
s(t) denotes the pulse envelope, and A weights the impor- 
tance of the penalty on the total field fluence. Solutions 
to the optimal control problem correspond to = 0. 
The problem of maximizing the expectation value of an 
observable corresponds to: 



$ X (C7) = Tr(U(T)p tf(T)Q) 



(2) 



where po is the initial density matrix of the system and 
G is an arbitrary observable operator, and the problem 
of maximizing the fidelity of a dynamical transformation 
W then corresponds to @|: 



$ 2 ([/) = Re TriAW^UiT)) 



(3) 



where W is the target unitary transformation and A is 
any Hermitian matrix. In the present study, we ignore 
the fluence penalty term because its effect is system- 
specific and not revealing with respect to the universal 
complexity of quantum control problems. 

We are interested in the convergence to the global 
optimum of the gradient flow trajectories of the objec- 
tive functions, which are the solutions to the differential 
equations 

^) =vmt))=a ^n (4 , 

where s is a continuous variable parametrizing the algo- 
rithmic time evolution of the search trajectory and a is 
an arbitrary scalar that we will set to 1. The complexity 



of quantum control optimization is associated with the 
scaling with system size of the rate of convergence of the 
function 



/ da' [ dt 


'$1,2(S',T)" 


loo JO 


Se{s',t) 



$i, 2 (s ) T)-$ li2 (oo,T) = a f da' f 

Joo Jo 

to zero. The gradients ^ can be shown to be, respec- 
tively us m G2: 



= Tr([e(T),poM*)), 



(5) 



and 



6<&2 

5e(t) 



Tr ((AAWU - U f WAA r )p(t)) , (6) 



within the electric dipole approximation, where /j,(t) = 
— j;U^(t, 0)liU (t, 0) is the time-evolved dipole operator of 
the quantum system and we have adopted the shorthand 
notation U=U(T). 

The gradient on e(t) is related to the gradient on U (N) 
through 



(5$ _ SUij d$ 



(7) 



Now suppose that we have the gradient flow of e(s,t) 
that follows ([U and let U(s), the system propagator at 
time T driven by e(s, t), be the projected trajectory on 
the unitary group hi (N) . The (algorithmic) time deriva- 
tive of U(s) is then 



dJ7y(s) 



T SUijis) de{s,t) 



ds J Se(s,t) ds 
which, combined with (01 and ([7]), gives 



dt 



ds 



tiUi 3 (s) ^ SU pq (s) d$ 



E 



o Se(s,t) Se(s,t) dU p 



(8) 



-dt. (9) 



It is convenient to write this equation in vector form, re- 
placing the N x N matrix U(s) with the N 2 dimensional 
vector u(s): 



du(s) 
ds 



<5u(s) Su T {s) 

Se(s,t) Se(s,t) 



dt 



V$ [u(s) 



G[e(s,t)]V$[u(s)] (10) 



where the superscript T denotes the transpose. This 
relation implies that the variation of the propagator in 
U(N) caused by the natural gradient flow in the space 
of control fields is Hamiltonian-dependent, where the 
influence of the Hamiltonian is contained in the 7V 2 - 
dimensional symmetric matrix G[e(s,t)]. Within the 
electric dipole approximation, G is given by 



Gij, pq {a) = I pij{s,t)p pq {s,t)dt. 
o 
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As we will discuss below, it is possible to eliminate 
G in this differential equation, and hence to explicitly 
follow the U(T) gradient flow, by adopting an alterna- 
tive algorithmic step on e(t). Importantly, the prop- 
erties of the map e{t) — ► U(T) for finite-dimensional 
quantum systems render this fairly simple to achieve 
across many diverse systems, albeit with an error that 
is Hamiltonian-dependent. For the present purpose, we 
can effectively ignore G because we are interested in 
comparing the complexity scalings of the flow trajec- 
tories for gate and observable control problems on the 
same quantum system. Averaged over many initial con- 
ditions, the time required for convergence of the U(T) 
gradient flows will determine this comparative scaling. 

Because the kinematic gradient flows associated with 
these differential equations evolve on the continuous 
space of a Lie group (for $2) or its adjoint orbit of skew- 
Hermitian matrices (for $1), the optimization processes 
can be considered to be analog rather than discrete com- 
putations. In prior work [12j . we demonstrated that 
these flow trajectories can in fact be viewed as dynam- 
ical systems whose equations of motion guide the al- 
gorithmic time evolution of gradient-based algorithms. 
We showed that although the infinite-dimensional gra- 
dient flow equations do not admit analytical solutions, 
the equations of motion for the gradient flow lines on 
the unitary group domain, namely 



HZ) =A<-UAU 
ds ) 2 



can in principle be exactly solved and explicitly inte- 
grated them for specific cases of interest. The solutions 
to the optimal control problem are the critical points of 
the flow, i.e. V$(£7) = 0, or equivalently, the equilibria 
of the dynamical system. In this work, we use the inte- 
grated gradient flow equations to derive upper bounds on 
the time for convergence to the solutions of these prob- 
lems and thereby assign the problems to computational 
complexity classes. 

The possible outputs of the dynamical system exe- 
cuting the analog computation are the attracting fixed 
points of the gradient flow of the control objective func- 
tional. Just as a discrete computation can be assigned 
a measure of time complexity, i.e. the scaling with sys- 
tem size (e.g., polynomial or exponential) of the time 
required to solve the problem, analog computations can 
be assigned time complexities, which correspond to the 
scaling with system size of the time required for the dy- 
namical system to converge to a vicinity of the optimal 
fixed point. Time complexity of an analog computation 
can be analytically determined for systems that are com- 
pletely integrable; i.e., for those systems whose equa- 
tions of motion can be exactly solved. For this purpose, 
we apply the recently developed theory of complexity for 
continuous time dynamics (T3 , ITij . 



For $1, the gradient flow evolves on a polytope whose 
dimension varies depending on the spectrum of po- Since 
the flow is cubic in U, many trajectories U±(s) corre- 
spond to the same function J(s). Rewriting the gradient 
flow as a quadratic function on the domain of Hermitian 
operators pr = U(T)p JJ Ji {T) pjj], one can solve explic- 
itly for a unique solution pt(s) corresponding to J(s). In 
the case that po has only one nonzero eigenvalue, corre- 
sponding to an initial pure state, we showed that under 
the change of variables pr{s) = \ip(s)){ip(s)), \ip(s)) — 
(ci(s),--- ,c N (a)), x(s) ee (|ci(s)| 2 ,- • • ,\c N {s)\ 2 "l_ the 
gradient flow can be explicitly integrated to give [12 ]: 



x{s) 



(| Ci (o)|V--,mo)| 2 ) 



ELM0)| 2 e 

ci(0)|V 



2s\i 



„2s\ r 



\cn(0)\ 2 



2 B 2s\i 



(11) 

(12) 



where Ai, • • • , Xn denote the eigenvalues of 9. The ex- 
plicit solution for the gradient trajectory of objective 
functional $2 was shown to be 

W^U(s) = (sinh(As) + cosh(As)VF t [/ )- 

( cosh(As) + smh(As)W^Uo)- 1 (13) 

where the initial condition is Uq = U(0) [l2j]. 

In the case of objective functional $1, the input for 
the computation may be viewed as being the eigenvalues 
of the observable operator, while the initial density ma- 
trix determines the initial condition. For the purpose of 
assigning a complexity class to the problem, the latter is 
generally taken as fixed. The time complexity is then ex- 
pressed as a function of the eigenvalues of 8. However, 
as discussed below, the complexity of the optimal con- 
trol problem can vary sharply as a function of the initial 
state as well. For objective functional $2, we consider 
the input for the computation to be the target unitary 
transformation W. The time complexity is expressed in 
terms of the eigenvalues of W. Again, time complexity 
can be modulated by varying the initial condition of the 
search away from the identity matrix I. 

It is possible to assess the convergence of the opti- 
mal control search in terms of the difference in objec- 
tive function values 3>(s) — $(00), but because of the 
degeneracy of solutions corresponding to a given value 
of 4>, a more precise assessment of complexity can be 
achieved in terms of the convergence of the distance be- 
tween the current guess and the global solution to the 
problem. As mentioned, the distance | \E(s) — E(oo)\ | 2 = 
Jq E(s,t)E(oo,t)dt is not an appropriate measure be- 
cause it is highly system-specific (precluding uniform 
behavior as a function of system size) and obscures the 
underlying geometry of convergence since $ is not ex- 
plicitly a function of E(t). The most appropriate choice 
is a distance on the space of solutions that displays the 
lowest degeneracy for a given value of $, while remain- 
ing independent of the system Hamiltonian. For <f>i, 



4 



this distance is ||f/(s) — f/(oo)||p, whereas for <j> 2 , it is 
the Euclidean distance on the polytope wherein the flow 
evolves. 

The attracting region R of an attracting fixed point 
is the subset of phase space wherein the distance to 
the point is monotonically decreasing with time. We 
consider the quantum control optimization problem to 
be solved to a desired precision, and the computation 
halted, if the gradient flow trajectory enters within an 
Cp-radius of the global optimum and is also within its 
attracting region. This definition ensures that for in- 
puts on which the e p -vicinity of the attractor is larger 
than the attracting region, entrance into the attracting 
region is required for halting of the computation. The 
computational complexity of the control optimization al- 
gorithm is then defined as the scaling with system size 
of the convergence time t c (H) — max [t c (e), t c (R)] to 
the intersection of the attracting region and the ball of 
radius e p , starting from a given initial condition. 



III. ANALOG COMPLEXITY OF QUANTUM 
CONTROL OPTIMIZATION 

A. Quantum observable maximization 

For optimization of $1 starting from an initial pure 
state, we can establish a bound on t c {e) as follows. We 
assume without loss of generality that the eigenvalues of 
O are arranged such that Ai>A2>A3>--->Ajv. For 
nondegenerate 0, the global optimum is then the unit 
vector ei» = e\ and the distance of the current point 
x(s) on the search trajectory to the global optimum can 
be written: 



| \x(s) — x(oo) 



IK«)H S 



(e 2 ^|ci(0)lV-- ,e 2 ^| Cjv (0)| 2 )- ei | i 



<2-2 r Al ' C ^' 2 (14) 



Eti^|c,(o)p 



since J2iLi x i = Ei=i c ? — 1 an d therefore YliLi x 



N ^2 



J2iLi c t — 1- Defining [i = Ai — A 2 , we then have 



jv 



^e 2s(A *- Al) | Cl (0)| 2 < e- 2 ^ + |ci(0)| : 



bound on the convergence time to the global optimum: 
ln(^«) 



tc(e) < 



2/t 



1 { 2N\ 



for e small. Note that the time scale for convergence to 
the global optimum (the linearization of the coefficient in 
the exponential) is given by fi. This is equal to the low- 
est eigenvalue of the Hessian of the objective function, 
H A (p(oo)) = - E 3 Yi( a ?i + ^ji)( A J~ Al )' where a ifc ,/3 ifc 
are the real, complex parts of an arbitrary Hermitian 
matrix A, obtained in [8( 1 . If the maximum eigenvalue 
of is degenerate with multiplicity fc, such that Ai = 
A 2 = • • • = Afe, and A& > Aj, j = k + 1, • • • , N, then 
the dynamics converges to the point ^(1, • • • , 1, 0, ■ • • ,0) 
[l2| . In this case, the distance to the global optimum be- 
comes 

|| a =(»)-i(i,---.i,o,---,o)|| 2 < 

.-^^ MP , (15) 



2^=1 



corresponding to 



1 , (2Nk\ 

A bound on the convergence time to the attracting re- 
gion of the solution, t c (R), can be derived as follows for 
the general case of a fc-fold degenerate maximal eigen- 
value of O. Again, without loss of generality, we assume 
Ai^..^ = A(!) > Afe + i > • • • > Xn- In the present case, 
t c (R) — t c such that for all s > t c , 

d 

— \\x{s) -x(oo)|| 2 = 
k N 

e 2sXl £ M0)| 2 ]T [(A (1) - A.Oe^lc^O)! 2 ] < 0. 

i=l j=l 

(16) 

This corresponds to the condition xi < for i > k 
[l2T ]. This condition holds iff Ej=i ^j x j > Aj, i — 
k + 1 , • • • , N. Insertion of the analytical solution above 
gives 

JV N 

£A,| Cj (0)|V^>A fc+1 £| Cj (0)|V^. 

3=1 3=1 

For the purposes of obtaining a bound on t c , we can 
rewrite this 



such that the bound on the distance to the solution be- 
comes 

\\x(s) ~ x(oo)\\ 2 < 2 - 2(1 + e- 2 ^ s |c 1 (0)|- 2 )- 1 . 

For simplicity, we choose the identity vector -4(1, 1) 
as the initial state. More generally, Xi*(0) will scale in- 
versely with N for randomly chosen initial conditions on 
the interior of the simplex. We then obtain the following 



£(A (1) |c 2 (0)| 2 -A fc+1 | Cfc+1 (0)| 2 ) e A <^= 



> 



(N-k- 2)A fe+1 |c fe+1 (0)|V fc + 1 \ (17) 



1 It can be shown that the expression derived in Q for "HaC^oo)) 
is equivalent to that for "H a(p{oo)) . 
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which can be solved for t c (R) to give: 
1 



t c ,i(R) < 



In. 



(1) 



A 



fc+i 



(N -k-2)\ k+1 \c k+1 (0)\< 



Eti [A(i ) |c l (0)P-A fe+1 | Cfc+1 (0)|^ 



(18) 



Therefore, the attracting region of the solution x(oo) = 
• • ■ , 1, 0, • • • ,0) does not cover the entire domain. 
Again taking the initial condition ^(1, 1), this bound 
becomes t cA (R) < 1 In ^zg^y , The upper bound 
on the computation time of the problem is then given 
by 



i c ,i(f0 = max [i c ,i(e)i*c,i(-R) ] < 



1 

2/x 



, /2Nk\ 
In - 



) 



2 In 



(N — k — 2)A fc+1 



fc(A(i) - Afe+i) 



(19) 



B. Quantum gate control 

For objective function $2, the control optimization 
is generally initiated with Uq at or near the identity 
transformation In, such that W'Uq = W'. Writing 
U'(s) = W^U(s), the solution to the problem then cor- 
responds to U'(oo) = In. We consider the special case 
where the Hermitian matrix A = I: 

U'(s) = (sinh(s/)+cosh(s/)W^ t )(cosh(s/)+sinh(s/)iy t )" 1 

For $2, the natural distance measure for assessment of 
convergence is given by the Frobenius matrix norm be- 
tween the current transformation U'(s) and the solution 
In: 

\\U'(s) - [/'(oo)||! = Tr[(U'(s) - I N )HU'( S ) - I N )} 

Diagonalizing W via the unitary transformation V, we 
have 

W = V^AV 



where A is the matrix of eigenvalues of W, 



A = 



/ 



(20) 



We previously showed [12j that this transformation al- 
lows the expression above to be simplified as 



N 



\\U'(s) 



IN 



k=l 



tanhs + e~ i0k 
~ lSk tanh s 



- 1 



N 

E 



2(1 — cos 9k) , , ,„ , 

i ^ g-C 1 - tanhs) 2 (21) 

1 + 2 tanh s cos 6 k + tanh s 



The partial derivative of each term in this expression 
with respect to cos 9 k monotonically decreases with 
cos 9k , so the greatest term in the sum will correspond 
to the lowest value of cos#fc, denoted cos#o- Therefore, 
the sum is bounded from above by A-times the term 
corresponding to 9q, 



\\U'( S )-I N \\ 2 <N 



2(1 - cos6> ) 



1 + 2 tanh s cos 9q + tanh s 



(1-tanhs) 2 



Under the change of variables x = tanhs, the condition 
for entering an e-vicinity of the solution is 



A 



2(1- cos g ) f1 r]2 
1 + 2x cos 9q + x 2 



-{l-xf <e. 



To derive tight upper bounds on x c and hence t c , 
we must therefore solve the quadratic equation (2A — 
2 A cos 9 - e)xl max + (4 A cos 9 - 4 A - 2e cos 9 )x c , max + 
(2A — 2 A cos 9q — e) = 0. The solutions to this equation 
are: 



= tanh t c 



4 A - 4 A cos 6»o + 2e cos 9 ± y/ (e 2 - 4 Ae) (cos 2 9 - 1 ) 



4A - 4A cos O - 2e 
I 



One solution gives identically a+ maX:+ = 1, correspond- 
ing to i c ,max,+ = 00, which only provides an upper 
bound on the convergence time. The other solution gives 
%cmax,- < 1, *c,max- < 00. In order to study the 



scaling of t c 



with the dimension of the system, we 



5 = 1 - x 



c.max. — 



2e(l - cos0 o ) + yj{e 2 - 4Ae)(cos 2 O - 1) ^ 
4A(l-cos0 o )+2e 

Sing ° 1 /I=o 1 /I (22) 
l-cos0oVA V A 



G 



for small e. This corresponds to 

t c {e) < t c ,ma X ,-(e) = In ( jj-jj a:c ' max '~ ^ = 

For $2, the distance to the solution U'(T) = In, 
U(T) — W* is a mqnotonically decreasing function of 
algorithmic time s [12], and hence in this case, the at- 
tracting region of the solution covers the entire manifold 
U(N) and t c (R) — 0. Thus, an upper bound on the com- 
putation time for the problem is 

1 /4N\ 

i c , 2 (H) = i c ,2(6)<-ln(— ) (24) 

for small e. 



C. Assignment to complexity classes 

A continuous time problem is said to be in the com- 
plexity class CLOG (continuous log) if it has a polyno- 
mial number of variables (here, the system dimension) 
and a logarithmic time complexity. Kinematic optimiza- 
tion of objective functions $i and <&2 therefore belong to 
the complexity class CLOG. CLOG is the analog coun- 
terpart of the classical complexity class NCi, the class 
of problems that can be efficiently solved on a discrete 
parallel computer, meaning problems that are decidable 
in polylogarithmic time on a parallel computer with a 
polynomial number of processors. It is the lowest time 
complexity class, lying immediately below P (polyno- 
mial time complexity). For real control optimization 
algorithms evolving in discrete time, the kinematic com- 
ponent of these problems will have NCi complexity. 

Quantum optimal control problems aimed at maxi- 
mizing the expectation value of any observable operator 
starting from a pure state of the system all belong to the 
same analog time complexity class. However, their char- 
acteristic time scale for exponential convergence, when 
the optimization algorithm follows the gradient flow tra- 
jectories as faithfully as possible, differs in a predictable 
fashion. Specifically, the characteristic time scale for 
exponential convergence is given by min^i \\j — Ai|, 
i.e., the magnitude of the difference between the largest 
and second largest eigenvalue of the observable opera- 
tor. This corresponds to the magnitude of the smallest 
Hessian eigenvalue of the objective functional near the 
solution. As such, the rate of convergence for the prob- 
lem of driving a pure initial state to a pure final state 

- which corresponds to Ai = 1, A; = Vi > 1 above 

- is the greatest, such that this problem has the lowest 
computational complexity of all observable maximiza- 
tion problems. The complexity classes corresponding to 
observable maximization problems starting from general 
mixed states may differ and are expected to be higher, 
as evidenced by 1) the factorial (vs exponential) scaling 



of the number of critical manifolds [8[ and 2) the fact 
that the grad ient flow evolves on a higher dimensional 
polytope [lll.ll2l|. 

Similarly, the kinematic component of the problem of 
optimizing controls for the implementation of a quantum 
gate belongs to the complexity class CLOG irrespective 
of the identity of the target gate W. In contrast to ob- 
servable expectation value maximization, however, the 
characteristic time scale for exponential convergence to 
the solution of quantum gate optimization does not vary 
directly as a function of the particular incarnation of the 
problem embodied in the choice of W, since the eigen- 
values of the Hessian matrix do not vary as a function 
of the eigenvalues of W @ . Nonetheless, the eigenvalue 
spectrum of W does affect the convergence time indi- 
rectly in conjunction with the initial guess Uq. 

Although the problems belong to the same complexity 
class with respect to the scaling of the number of itera- 
tions required as a function of system size, they display 
distinct behavior in other regards. In particular, it can 
be shown [l2j that for both problems, there exists a set 
of initial conditions from which the flow does not con- 
verge to the solution. For $2, this corresponds to the 
case where cos 8q = —1, resulting in x c - = 1, i c ,_ = oo. 
Thus, if W has at least one eigenvalue equal to —1, the 
optimization does not converge if U$ is chosen as /, since 
U'(0) — W^Uq = . For $i, initial conditions where 
c.;* = (i.e., where the initial state resides on the so- 
called basin boundary of the simplex) do not converge. 
The scaling of the convergence time as a function of the 
distance to this pathological initial condition differs for 
the two problems [l2l ]. In these cases, the initial guess 
Uq should be modulated to facilitate the convergence of 
gradient algorithms. 

The complexity of optimal control landscapes for the 
optimization of quantum gates scales linearly with the 
number of qubits. If the problem size is measured in 
qubits, optimization of $2 belongs to the complexity 
class CP (continuous P), the set of problems with a poly- 
nomial number of variables and polynomial time com- 
plexity on a continuous variable computer. CP is the 
analog counterpart of the discrete complexity class P. 
In previous work aimed at optimizing quantum gates 
through OCT, the scaling of the number of control field 
iterations was reported to be exponential in the number 
of qubits [HI]. In these studies, the gradient flow was not 
followed directly by the optimization algorithm. The re- 
sults here suggest that exponential speedup should in 
principle be possible for the implementation of quantum 
gates through the use of algorithms that follow the gra- 
dient flow on the unitary group. 

As such, a central issue concerns the accessibility of 
these lower bounds on scaling through the design of 
gradient-based algorithms in both the laboratory and in 
simulations. These continuous time complexities repre- 
sent lower bounds on the complexities actually achieved 
by discretized Euclidean search algorithms that do not 
precisely follow the Riemannian curvature of the gradi- 
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ent flow trajectories. The increase in complexity induced 
by (Euclidean) discretization depends on the curvature 
of the trajectory. A quantitative analysis is beyond the 
scope of the current work, but the path length of the gra- 
dient flow trajectory offers insight into its geometry. On 
the unitary group U (N) , the path length of the gradient 
flow trajectory is given by 

L(U ,t c )= r JTv(Ut( s )U(s))d Sl 
Jo 



where the Riemannian metric on the unitary group is 
defined as (A, B) — ^(A^B). For the optimization <f>2, 
this integral can be shown to be 

L(U , t c) = 2 J o C ^N-ReTr [{WW) 2 ]ds, 
which can be expressed as a sum of N scalars: 



L(U ,t c )=2 ds A 
Jo 



N 



2 sin 2 6» fc (l-tanh 2 s) 



\ (1 + tanh 2 s) 2 + 4 tanh s [cos 2 6 k + (1 + tanh 2 s) cos 9 k ] 



Again setting x = tanhi, the integrand can be further 
simplified using the relation 



2Nsin 2 fc (l -x 2 ) 



2\2 



(1 + x 2 ) 2 +Ax [cos 2 6> fe + (1 + x 2 )cos6'fc 



< 





(l + X 2 )y/T-X 
dx 



din 



(25) 



1+x 
1-x 



(1 + x 2 )Vl - x 
sin 6»| f Xc dx 



VI - %~c Jo 1 + x 2 
sm8\ arctana: c 



sin 01 



VI - x c 



sin 



0| e -V4 iV 3/4 > 



Therefore, we see that the path length - which unlike the 
convergence time, does not depend directly on the mag- 
nitude of the gradient and hence is less sensitive to the 
effects of discretization - scales approximately linearly 
with N, roughly consistent with the scaling observed in 
kinematic simulations of this objective [l6[. 



IV. DISCUSSION 

The complexity results derived herein are based on 
analysis of the kinematic contribution to optimal control 
gradient flows, and have neglected the dynamical contri- 
bution of the e(t) — ► U(T) map to the scaling of search 
effort. For an arbitrary Hamiltonian, if one averages over 



many initial conditions and target gates/states, the rel- 
ative convergence times of gate and state control opti- 
mization should correspond to those of the kinematic 
gradient flows |12j. However, since the Hamiltonian, 
and hence the properties of the e(t) — ■> U(T) map must 
change with system dimension, the absolute scaling for 
either problem may in principle deviate from that pre- 
dicted here. The primary importance of the above re- 
sults is that the relative complexity scaling of these two 
problems will nonetheless be identical, since the contri- 
butions of the functional derivatives to scaling, 
averaged over initial conditions and target gates, will be 
the same for each. As such, quantum pure state and 
gate control problems do belong to the same complexity 
class. 

From a practical perspective, since the kinematic com- 
plexity of these problems is as low as possible, the scal- 
ing of control search effort should not prohibit the effi- 
cient application of gradient-based algorithms to high- 
dimensional systems. For state-to-state coherence trans- 
fer problems, this prediction is borne out by both sim- 
ulations and experimental evidence [1 €31 ] - However , the 
reported scaling of gate optimization simulations 0] is 
worse than that predicted according to the present the- 
ory, suggesting that a reexamination of current quantum 
gate control algorithms is warranted. For control prob- 
lems whose kinematic complexity class exceeds CLOG, 
it may be advantageous to apply matrix tracking algo- 
rithms (below) rather than gradient-based algorithms 
in order to achieve the above lower bounds. Evidence 
suggests that observable maximization problems start- 
ing from an initial nondegenerate mixed state belong to 
this category [l7j . 

For general po and O, the gradient flow for objective 
function $1 is an " isospectral" flow [H, OJ, HI • In N di- 
mensions, this flow has N integrals of the motion that are 
in involution, which is the classical definition of complete 
integrability for a dynamical system. From the point of 
view of the modern theory of integrable systems, the 
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double bracket flow can be shown to represent a type of 
Lax pair, a general form that can be adopted by all com- 
pletely integrable dynamical systems 21]. As such, the 
system-independent complexity of all discrete quantum 
observable maximization problems can in principle be 
analytically determined [12]. In contrast, the kinematic 
gradient flows of the gate and observable control prob- 
lems in classical mechanics [3, [22| are not integrable, and 
hence cannot be assigned to analytic complexity classes. 

In the present work, we have focused attention on 
the complexity of optimizations that employ local gra- 
dient search algorithms. The favorable scaling of the 
convergence times for kinematic gradient flows, and the 
scarcity of abnormal extremals in discrete quantum con- 



trol problems [4j motivate a rigorous definition of the 
universal complexity of quantum control problems in 
terms of the scaling of the expense of tracking such 
kinematic paths through elimination of the Hamiltonian- 
dependent matrix G in the algorithmic step (TlT)|) . These 
global algorithms require an additional computational 
overhead of N 4 in order to invert the matrix G, as 
well as overhead for statistical estimation of the states 
or dynamical propagators based on experimental ob- 
servations 0]. Assignment of quantum control prob- 
lems to complexity classes based on the convergence of 
such global algorithms, which involves extensive sam- 
pling over Hamiltonian space, is the subject of a separate 
work [H[. 
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